% Plot ESR data
% hyzhou, 12/19/2017
% param stores values and corresponding errors

function [param,h] = plot_esr(freq, esr, fontsize)

h = plot(freq/1e9,esr);
hold on

fit_type = 'singleLor';

switch fit_type
    case 'doubleLor'
        % Fit with two Lorentzians
        fit_func = fittype('1 - a1./(1+((x-b+e)/c).^2) - a2./(1+((x-b-e)/d).^2)');
        [fit_res,~] = fit(freq'/1e9,esr', fit_func,...
            'startpoint', [(1-min(esr))/2 (1-min(esr))/2 2.87 0.01 0.01 0.001], 'lower', [0.00 0.00 2.4 0 0 0],...
            'upper', [1 1 2.98 0.1 0.1 0.1],  'robust', 'BiSquare');
        errs = diff(confint(fit_res),1)/4;
        plot(freq/1e9, fit_res(freq/1e9))
        xlabel('\omega (GHz)', 'fontsize', fontsize); set(gca, 'fontsize', fontsize);
        ylabel('Contrast', 'fontsize', fontsize)
        title(['Left: ',num2str(fit_res.a1*100,'%.2f'),'%, HWHM ',...
            num2str(1e3*fit_res.c,'%.1f'),' MHz, ',...
            'Right: ',num2str(fit_res.a2*100,'%.2f'),'%, HWHM ',...
            num2str(1e3*fit_res.d,'%.1f'),' MHz', sprintf('\n'),...
            'Center: ',num2str(fit_res.b*1000,'%.2f'),'\pm',num2str(errs(3)*1000,'%.2f'),...
            ' MHz, Half split: ',...
            num2str(fit_res.e,'%.5f'),' GHz'], 'fontsize', fontsize-4,...
            'interpreter','tex')
        
    case 'singleLor'
        % Fit with a single Lorentzian
        fit_func = fittype('1 - a./(1+((x-b)/c).^2)');
        [fit_res,~] = fit(freq'/1e9,esr', fit_func,...
            'startpoint', [(1-min(esr)) mean(freq/1e9) 0.01], 'lower', [0.00 2.4 0],...
            'upper', [1 2.98 0.1],  'robust', 'BiSquare');
        errs = diff(confint(fit_res),1)/4;
        plot(freq/1e9, fit_res(freq/1e9))
        xlabel('\omega (GHz)', 'fontsize', fontsize); set(gca, 'fontsize', fontsize);
        ylabel('Contrast', 'fontsize', fontsize)
        title(['Amp: ',num2str(fit_res.a*100,'%.2f'),'%, HWHM ',...
            num2str(1e3*fit_res.c,'%.1f'),' MHz, ',...
            'Center: ',num2str(fit_res.b*1000,'%.2f'),'\pm',num2str(errs(3)*1000,'%.2f'),...
            ' MHz'], 'fontsize', fontsize-4,...
            'interpreter','tex')
end
param = [coeffvalues(fit_res);errs];